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Abstract. This paper deals with the numerical reconstruction of the plasma current density 
in a Tokamak and of its equilibrium. The problem consists in the identification of a non-linear 
source in the 2D Grad-Shafranov equation, which governs the axisymmetric equilibrium of a 
plasma in a Tokamak. The experimental measurements that enable this identification are the 
magnetics on the vacuum vessel, but also polarimetric and interferometric measures on several 
chords, as well as motional Stark effect or pressure measurements. The reconstruction can be 
obtained in real-time using a finite element method, a non-linear fixed-point algorithm and a 
least-square optimization procedure. 



1. Introduction 

The problem of the equilibrium of a plasma in a Tokamak is a free boundary problem in which 
the plasma boundary is defined either by its contact with a limiter or as being a magnetic 
separatrix. Inside the plasma, the equilibrium equation in an axisymmetric configuration is 
called Grad-Shafranov equation [U El E]- The right-hand side of this equation is a non-linear 
source which represents the toroidal component of the plasma current density. 

An important problem is the identification of this non-linearity [H [6]. The aim of this 
paper is to present a method for real-time identification from experimental measurements, such 
as magnetic measurements on the vacuum vessel, polarimetric measurements (integrals of the 
magnetic field over several chords), MSE (Motional Stark Effect) and pressure measurements. 
The pressure is supposed to be isotropic. For the anisotropic pressure case, one can refer to [7|. 

The next section is devoted to the mathematical modelling of the equilibrium problem in 
axisymmetric configurations. The inverse reconstruction problem is adressed in the last section. 



2. Mathematical modelling of axisymmetric equilibrium of the plasma in a 
Tokamak 

The equations which govern the equilibrium of a plasma in the presence of a magnetic field are 
on the one hand Maxwell's equations and on the other hand the equilibrium equations for the 
plasma itself. 

The magnetostatic Maxwell's equations as follows are satisfied in the whole of space (including 
the plasma): 
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where B represents the magnetic field, u is the magnetic permeability and j is the current 
density. The first relation of (pQ) is the equation of conservation of magnetic induction and the 
second one is Ampere's Theorem. 

The momentum equation for a plasma is 



P ~dt = J X P 



(2) 



where u represents the mean velocity of particles and p the mass density. At the resistive 
time-scale the first term can be neglected [S] and the equilibrium equation for the plasma is 



Vp = j x B 



(3) 



This equation ([3]) means that the plasma is in equilibrium when the force Vp due the kinetic 
pressure p is equal to the Lorentz force of the magnetic pressure j x B. We deduce immediately 
from Q that 

B ■ Vp = (4) 



j ■ Vp = 



(5) 



Thus for a plasma in equilibrium the field lines and the current lines lie on isobaric surfaces 
(p = const.); these surfaces, generated by the field lines, are called magnetic surfaces. In order 
for them to remain within a bounded volume of space it is necessary that they have a toroidal 
topology. These surfaces form a family of nested tori. The innermost torus degenerates into a 
curve which is called the magnetic axis. 

In a cylindrical coordinate system (r, z, eft) (where r = is the major axis of the torus) the 
hypothesis of axial symmetry consists in assuming that the magnetic field B is independent of the 
toroidal angle (f). The magnetic field can be decomposed as B = B p + Ba, where B p = (B r , B z ) 
is the poloidal component and B^ is the toroidal component. From equation ([T]) one can define 
the poloidal flux vp(r, z) such that 
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Concerning the toroidal component we define / by 



B^ = J -t 



(7) 



where e<f, is the unit vector in the toroidal direction, and / is the diamagnetic function, 
magnetic field can be written as: 
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According to ([8]), in an axisymmetric configuration the magnetic surfaces are generated by the 
rotation of the flux lines ip = const, around the axis r = of the torus. 

From ([8]) and the second relation of (pQ) we obtain the following expression for j: 

j = j P + H 
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where j p and are the poloidal and toroidal components respectively of j, and the operator 
A* is defined by 

did. did. 
A*. = —(——) + —(— — ) (10) 
dr \xr dr dz fir dz 

Expressions (JS|) and ([!]) for B and j are valid in the whole of space since they involve only 
Maxwell's equations and the hypothesis of axisymmetry. 

In the plasma region, relation @ implies that Vp and are colinear, and therefore p is 
constant on each magnetic surface. This can be denoted by 

p = pW (ii) 

Relation ([5]) combined with the expression ([9j) implies that V/ and Vp are colinear, and therefore 
/ is likewise constant on each magnetic surface 

/ = m (12) 

The equilibrium relation ([3]) combined with the expression ([5]) and © for B and j implies that: 

Vp = -^V^ - -L^Vf (13) 

which leads to the so-called Grad-Shafranov equilibrium equation: 

-A*^ = rp'^) + — (//0W0 (14) 

where A* is the linear elliptic operator given by ()10|) in which fi is equal to the magnetic 
permeability of the vacuum. 

From ([9]) it is clear that right-hand side of (|14H represents the toroidal component of the 
plasma current density. It involves functions p{ip) and f(ijj) which are not directly measured 
inside the plasma. 

In the vacuum, the magnetic flux tp satisfies 

-A> = (15) 

The equilibrium of a plasma in a domain $7 representing the vacuum region is a free boundary 
problem. The plasma free boundary is defined either by its contact with a limiter D (outermost 
flux line inside the limiter) or as being a magnetic separatrix (hyperbolic line with an X-point, 
X). The region Q p C 0, containing the plasma is defined as 



(16) 



Figure 1. Definition of the plasma boundary (thick blue line). Left, JET (Joint European 
Torus) example, X-point configuration. Right, TORE SUPRA (the CEA-EURATOM Tokamak 
at Cadarache) example, limiter configuration (the limiter is represented by the black line). The 
thin blue lines represent iso-contours of ip. 



where either ifo = maxtp in the limiter configuration or ip^ = ip(X) in the X-point configuration 
(see Fig. HJ 

Assuming Dirichlet boundary conditions, h, are given on T = dO, which is the poloidal cross- 
section of the vacuum vessel, the final equations governing the behaviour of ip(r,z) inside the 
vacuum vessel, are: 

-A*p = [rA$) + ifl$)]xnp in Q 

tp = h on r 

-t ip — max ip 

with AN)) = p'{4>) and B(ip) = —(ff')($), $ = G [o, 1] in Q p (this normalized flux 

Hq ipb — max ip 

is introduced so that A and B are defined on the fixed interval [0, 1]), xn p is the characteristic 
function of Cl p . 

The aim of the following section of this paper is to provide a method for the real-time 
identification of the plasma current i.e. the non-linear functions A and B in the elliptic equation 
(EFT 



3. The inverse problem 

3.1. Experimental measurements 

The given experimental measurements are: 

• the magnetic measurements 

— ip(Mi) = hi on r, given by the flux loops (see Fig. [2|). Thanks to an interpolation 
between points Mj these measurements provide the Dirichlet boundary condition h. 



1 dtp 

— ———(Ni) = gi on r, which corresponds to the component of the magnetic poloidal 
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field, measured by the magnetic probes (see Fig. [21), which is tangent to the vacuum 

vessel. Indeed from Eq. ([8]) the tangential component of B p is equal to the normal 

1^ ,1, 
component — - — oi —vip. 
r on r 

the polarimetric measurements which give the Faraday rotation of the angle of infrared 
radiation crossing the section of the plasma along several chords Cf. 

n e (4>) dip 
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where n e represents the electronic density which is approximately constant on each flux 

d 

line, B\\ is the component of the poloidal field tangent to C{ and — represents the normal 

1 On 
derivative of if) with respect to C,. 

• the interferometric measurements which give the density integrals over the chords Cj 

n e (4>)dl = (5i 

Ci 

the kinetic pressure measurements obtained from density and temperature measurements, 
for instance in the equatorial plane: 

p(r,0) =p d {r) 

and MSE (Motional Stark Effect) angle measurements taken at different points Xi = (rj, z\): 

m(B r ,B z ,B (j} ) i = ji 

with 

+ ( en d u \\ a\B r + a 2 B z + 
tan(m(B r ,B z ,B^)) 
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3.2. Statement of the inverse problem 

The numerical identification problem is formulated as a least-square minimization with a 
Tikhonov regularization. The cost function is defined as: 



with 



J(A, B, n e ) = J + K X J X + K 2 J 2 + K 3 J 3 + K 4 J 4 + J e (18) 
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Figure 2. Left: the straight green lines represents the chords used for polarimetry and 
interferometry measurements. Right: part of the vacuum vessel. At the bottom middle an 
example of finite element mesh used for numerical simulations (see next Section). 

and Ki, K2, -K3 and K4 are weighting parameters enabling to give more or less importance 
to the corresponding experimental measurements [5]. 

The inverse problem of the determination of A and B is ill-posed. Hence a regularization 
procedure can be used to transform it into a well-posed one [9]. The Tikhonov regularization 
term J e constrains the function A, B and n e to be smooth enough and reads: 

J e = ei [ 1 [A"(x)] 2 dx + e 2 ^ [B" (x)] 2 dx + e 3 {\nl(x)} 2 dx 
Jo Jo Jo 

where e±, 62 and £3 are the regularizing parameters. 

It should be noticed that the electronic density n e does not intervene in Eq. (|17p . However as 
soon as we want to use the polarimetric measurements it is necessary to include n e (and hence 
interferometry) in the identification procedure. The inverse problem can finally be formulated 

as, 

Find A*, B*, n* such that : 

(19) 

J{A*,B*,n* e ) = inf J{A,B,n e ) 
3.3. Numerical identification 

Problem (|17|) is solved using a finite element method [TP]. Let i/ 1 (0) and V = Hq(£2) denote 
the usual Sobolev spaces. The finite element approximation is based on the following weak 



formulation: 

Find if) € H 1 ^), such that ip = h on T, and 



1 /-,_!_ (20) 



Classically is approximated using triangles by a polygonal domain Q^, the space V is 
approximated by a space of finite dimension n. A PI finite element method is used, in 
which functions of Vh are affine over each triangle and continuous on the whole domain. 

Let K denote the finite element stiffness matrix. Let us also (abusively) denote by ifi G R n 
the components of the magnetic flux fonction approximated in V/j. 

The unknown functions A, B and n e are approximated by a decomposition in a reduced basis 

{4>i)i=l,...m 

A(x) = ^ai4>i(x) 

i 

B(x) = Y,bMx) 

i 

n e (x) = ^2ci4>i(x) 

i 

This basis can be made of different types of functions (polynomials, splines, wavelets, etc . . . ) 
[6]. Let u be the vector of R 3m defined by u = (a±, . . . , a m , bi,..., b m , c\, . . . , c m ). With these 
notations the discretization of problem (|20p can be written as follows: 



Given u € R 3m , solve the fixed — point equation 

(21) 

Ktp = D(ip)u + h 

Where D{ip) denotes the n x 3m "plasma current matrix", and K is the stiffness matrix 
modified in order to impose the Dirichlet boundary condition represented by h. 
The discrete inverse optimization problem is: 

Find u minimizing : 

J(u) = ||C(^-jfe|| 2 +« T A« (22) 



with ij) satisfying (I2ip 

where C(ip) is the observation operator. The quantity C(ip)ip represents the outputs of 
the model corresponding to the experimental measurements, given in a previous subsection, 
denoted by k. The matrix A represents the regularization terms. The first term of J in Eq. (|22p 
corresponds to Jo + K\ J\ + K2J2 + K3J3 + J4 and the second to J e . 

In order to solve this problem we use an iterative algorithm based on fixed-point iterations 
for Eq. (j2T|) and the normal equation of Eq. (|22l) . 



3.3.1. Algorithm At the re-th iteration, tp n and u n are given. The non- linear mapping u 1— » tp(u) 
is approximated by the affine relation 

i> = K- 1 [D{i> n )u + h] 

and the cost function to be minimized by 

J(u) = ||C(^ n )V> -k\\ 2 + u T Au 

= \\Caj n )k~ l D^ n )u + {Ci^K^h -k)\\ 2 + u T Au 

= \\E n u + F n \\ 2 + u T Au 



with obvious notations. The normal equation 



{E T n E n + A)u = -E T n F n 
is solved to update u n to u n+ \. Then a fixed-point iteration for Eq. (f2T1) enables the update of 

tp n to Vn+1 

tp n+1 = K^lD^lfj^Un+i + h]. 

Since the algorithm is usually initialized with the equilibrium at a previous time step, two or 
three fixed-point iterations are usually enough to ensure convergence. 

3.3.2. Equinox software Based on the algorithm presented above, a C++ software, called 
EQUINOX |114 [12l [T3] has been developed in collaboration with the Fusion Department at 
Cadarache, and has been implemented for JET (Joint European Torus) and for TORE SUPRA 
(the CEA-EURATOM Tokamak at Cadarache). Figure [3] shows a graphical output of Equinox. 
With all these techniques it is possible to follow the quasi-static evolution of the plasma 
equilibrium, either in TORE SUPRA or JET configurations, with free boundaries defined either 
by limiter contact or with an X-point. It is also possible to simulate ITER configurations. 




Figure 3. An output of EQUINOX. The plasma is in an X-point configuration. On the left 
column the identified p', ff and n e functions as well as the toroidal current density j and the 
safety factor q are displayed in terms of ip and of r (in the equatorial plane) . 



4. Conclusion 

We have presented an algorithm for the identification of the current density profile in Grad- 
Shafranov equation from experimental measurements. 

The decomposition of the unknown functions p'(ip) and ff'(ifi) in a reduced basis makes it 



possible to do the reconstruction in real-time. 

The choice of this reduced basis must still be improved and optimized (robustness, precision, 
...). 

Real-time reconstruction makes possible future real-time control of the current profile |14j . 
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